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0\ ■ Abstract 

a^ 

as 

The Debye- Waller factor and the mean-squared displacement from lattice 

O I sites for solid He and He were calculated with Path Integral Monte Carlo 

\m/ ' at temperatures between 5 K and 35 K, and densities between 38 nm~^ and 

^^ ■ 67 nm^'^. It was found that the mean-squared displacement exhibits finite-size 

^ . scaling consistent with a crossover between the quantum and classical limits 

of N^'^'^ and N~^'^, respectively. The temperature dependence appears to 

^ [ he T^, different than expected from harmonic theory. An anisotropic /c^ 

(^ ' term was also observed in the Debye- Waller factor, indicating the presence 

^J^ ■ of non-Gaussian corrections to the density distribution around lattice sites. 

Our results, extrapolated to the thermodynamic limit, agree well with recent 

On , values from scattering experiments. 

0\ 
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I. INTRODUCTION 



Solid helium at low temperatures is the prototype quantum crystallil, since its ground 



O ■ state is not well-described by harmonic perturbations about a minimum potential configura- 

^ ' tion. For temperatures between 5 K and 35 K, and densities between 38 nm~^ and 67 nm~^, 

the kinetic energy is always greater than about 25 K and dominates over thermal energies. 
/\^ • Solid helium has a wide range of experimentally accessible temperatures and densities. With 

j^ ■ new scattering sources, the density distribution and other correlation functions can be mea- 

sured with unprecedented accuracy, allowing for a careful comparison between theory and 
experiment. 

It is possible to calculate the properties of helium very accurately using Monte Carlo 
methods, because the inter-atomic potential is accurately known, more accurately than any 
other atomic or molecular condensed matter system. Additionally, for bosonic and distin- 
guishable particle systems. Path Integral Monte Carlo methods can calculate equilibrium 
properties directly from an assumed Hamiltonian without significant approximation. Even 
in the case of solid ^He, effects of Fermi statistics can be neglected for temperatures above 
0.1 K and densities slightly away from melting. 

In scattering, the Debye- Waller factor is the fractional intensity shift due to recoilless 
processes, and can be directly related to the one-dimensional mean-squared displacement 
of particles from their lattice sites (m^). In this paper we compute (u^) and compare it to 
experimental results, obtaining good agreement. In the process we make several interesting 
observations. First, the density distribution of helium atoms is slightly non-Gaussian, a 



fact that had been speculated on but not yet observed. Second, we observe some unusual 
dependence of {v?) on the size of the finite system being simulated, indicating a scaling 
crossover between the quantum and classical regimes. The dependence of {v?) on the number 
of atoms A^ being simulated appears to be well- described by harmonic theory. Finally, the 
temperature dependence appears to be T^, rather than the T^ predicted by harmonic theory. 



A. The Debye- Waller Factor 

The static structure factor, as measured in scattering, is defined as 

S{k) = (Ip.l^) (1) 

where pk = -^ S^*'''"'. The structure factor can be computed directly with PIMC for 

i 

finite systems of A^ ^ 10000 atoms. In a solid, the structure factor has large peaks at the 
reciprocal lattice vectors of the perfect lattice. It is usually assumed that the magnitude of 
S{k) behaves as 

S{k) oc exp(-2P^) ^ exp(-A;2(M2)) (2) 

where exp(— 2Vr) is the Debye-Waller factor. This equation relies on the assumption that 
the particle densities are normally distributed about the lattice sites. 

We will now derive a more general form of Eq. (^. We will assume a finite system, with 
periodic boundary conditions and precisely A^ particles at positions r^ and A^ perfect lattice 
sites Zj. Using particle symmetry, we can rewrite Eq. (0) as: 

5(A;) = l + (A^-l)(e''^-('-^-'-j)),-, (3) 

where the angle brackets ()j denote an average over both the thermal density matrix and 
particles j 7^ 1. When k is a reciprocal lattice vector, 

where Mj = k ■ (ri — Zi) is the displacement of particle i from its lattice site in the direction 
of k. We assume a simple Bravais lattice, and assign each particle to a lattice site. Now 
consider the variable x = Ui — Uj. Using the cumulant expansionB to evaluate the average 
of an exponential in terms of the moments of x, we can write: 

(e^^-)=exp(-^(x^) + ^((x^)-3(a:y) +...), (5) 

since the odd powers of x in the expansion vanish under the interchange 1 <-*> j allowed 
by particle symmetry. For a system much larger than the correlation length ^ of m, which 
is finite in solid helium, ui is uncorrelated with Uj except for the neighbors of particle 1. 
Hence: 

{x') = 2{u') - 2{mu,) - 2{u') + O {{i/Lf) 

{x^) = 2{u^) + Q{uY - A{u\u,) - A{uiu]) = 2{u^) + Q{uY + O {(^/Lf) . (6) 



Here, L is the box length. Combining Eqs. (^5,6), we obtain 

S{k) = 1 + (iV - l)exp {-k'^u") + ^{uY^ , (7) 

where the kurtosis k is defined as the relative deviation of the fourth moment from a normal 
distribution, 

It vanishes if the density distribution is normal in the scattering direction. Previous analysis 
of this kind has assumed Gaussian fiuctuations, and hence neglected all higher order terms 
beyond the first. However, we find that the kurtosis is not precisely zero and has directional 
dependence. 



II. PIMC CALCULATIONS OF SOLID HELIUM 

Path Integral Monte Carlo simulations were performed as discussed in ref.u. The system 
being simulated consists of particles (Boltzmannons) in a box with periodic boundary con- 
ditions at a fixed density. The helium atoms were assumed to interact pair-wise with the 
Aziz potentialu. Although small errors in energy are expected with this potential due to the 
absence of three-body interactions, we expect the pair potential to describe well the density 
distribution due to the fact that three-body and higher order contributions are smoothly 
varying with respect to atomic positions. We implicitly test this assumption by comparing 
to experimental values. The pair potential was set to zero for inter-atomic distances greater 
than 6 A. We determined that this is the minimum cut-off that could be used without caus- 
ing systematic errors in the mean-squared displacement. The pair potential was used to 
compute the exact pair density matrix for the system with an imaginary time-step equal to 
r = 1/160 K~^. Using this action, the time-step error was found to be negligible. We used 
neighbor lists to achieve linear scaling of the computer time versus the number of atoms, 
and were able to simulate systems of up to 3000 atoms. 

We computed the Debye- Waller factor two different ways. First, {v?) was computed 
directly from the distance of the atoms from their lattice sites, 

where M is the number of imaginary time slices, A^ is the number of particles, and Tjj is the 
position of particle i at imaginary time slice j. The factor of 1/3 arises because in a cubic 
crystal we can also average over the three spatial directions. 

For the sake of convenience, it is useful to forbid particle exchanges between lattice sites 
so that the particles do not have to be periodically re-assigned to the nearest site. We 
assured localization to a lattice site by "tethering" each particle: specifying a distance from 
its lattice site past which all attempted moves were rejected. The tether distance, 2.6 A, was 
chosen to be on the order of, but slightly less than, the average nearest-neighbor distance 
and did not introduce a noticeable change in (u^) or the structure factor. 



The second method for computing the Debye- Waller factor was to calculate S{k) directly 
from Eq. (|l|) from the set of all reciprocal lattice vectors oi k < 9 A~^ and then use Eq. (^ 
to determine (u^) and the kurtosis k by a least-squares fit to \og{{S{k) — 1)/{N — 1)). This 
method has the advantage of not requiring tethering or indeed any of the assumptions used 
in deriving Eq. (0). It also allows one to determine not only (u^) but also any non-Gaussian 
components. Agreement between the two approaches shows that correlations in u in Eq. (6) 
do decay rapidly. Calculating (u^) either directly with Eq. (^ or from fitting to S{k) always 
gave the same value within statistical error, with similar error bars (see Fig. [0]). 

III. FINITE-SIZE SCALING 

Before we can present the comparison to experiment, we must extrapolate the results 
obtained for a finite system to the thermodynamic limit. We observed a very slow conver- 
gence. To carefully examine the finite-size effects we simulated much larger systems (up to 
3000 atoms) than had been done previously with PIMC Even with the increased range of 
system size, the finite-size effects are not well- described by a power law, but are instead in a 
crossover region. For the temperature and density values studied in this paper, the crossover 
region appears to span the range of system sizes available to present computer simulation 
(A^crossover < 10^), requiring careful fitting to obtain values in the thermodynamic limit. 

Young and AldeiO have used Debye theory to analyze the finite-size dependence of a 
system of classical hard spheres and determined that (u^) oc p~'^^^TN~^/^. They found that 
this scaling was accurately able to fit values obtained with molecular dynamics simulations. 
Recent classical Monte Carlo calculation on a Lennard- Jones model of an fee solidH also 
found the same dependence. 

Runge and ChesteiQ looked at the size effects of a hard sphere system using PIMC. Using 
the same Debye theory but now taking quantum effects into account, they estimated that 
the finite-size effects scale as A^~^/^ at zero temperature, and found a crossover from the 
quantum to the classical scaling (A^~^/'^) as the temperature was increased. 

We used harmonic theory to derive a reasonable functional form for the finite-size effects 
of (v?) in a crossover region and the width of that crossover region. In general, one can 
write the mean-squared displacement asB 

3 [ S'k f duj A{k,uj) 



lu') = a-" / --^ / — :r'^ (10) 

^ ' J (27r)3 i 27r e/3^ - 1 ^ ' 

where y4(k, cu) is the spectral function for the displacement-displacement Green's function, 
a? = 1/p is the volume of the unit cell, and /? = l/ksT. It is useful to define 



so that 



{u') = j^Jd'kA{k). (12) 



Let us assume that the effect of the periodic boundary conditions is to replace the integral 
in Eq. (|T0|) with a sum: 

6{u') ^ {u')^ - {u')^ = (1^ / ^'k A{k) - ^ E ^(k) kl (13) 

The finite spacing of k is a function of the system size, and is given by kc = 2ti(p/NY^^ for 
a cubic simulation cell. The main contribution to the finite-size error is due to the omission 
of the k = value. If we further assume that for small values of k that A(k) factors into 
an analytic function (smooth and continuous near k = 0) and a singular factor \k\~'^, then 
it can be shown that the dominant term in Eq. (|T3|) is: 

6{u^) oc kl^" oc m"\ (14) 

According to this theory, one must determine the exponent z/ of the singular part of the 
Green's function A{k), at fc = 0. 
For a harmonic lattice: 

Aarmonic(k,U;) = [6{uJ - Uk) - S{u + LOk)], (15) 

zmujk 

where uJk is the frequency of a phonon of wavevector k and m is the mass. In the limit of 
small k, the dispersion Uk is linear in k. Integrating over uj, we obtain 

Aik) = ^^cotH^ (16) 



Thus for ksT ^ /icufc we find i^ = 2, but at T ^ 0, z/ = 1, as did Youngp and RungeQ. A 
crossover between the two scaling forms occurs at huk ~ ksT. 

Solid helium is known to have a significant phonon linewidth. For small k, the phonon 
linewidth can be approximated to be 7fc = 0.2co'fc^ii2l. For a damped harmonic lattice, 

Adan.ped(fc, '^) = ^^^^ ( (o; - c ^^ + 7^/4 " (c + c^.^ + 7^/4 j • ^^^^ 

Integrating over u, we find that for T = 0, 

A{k) = tan-M — ^ oc — , (18) 

mujk \ 7fc / UJk 



which gives v = 1. For T > 0, Eq. ([1^) can be integrated numerically. As before, we find 
z/ = 2 when ksT 3> hjJk and z/ = 1 when ksT -C hjJk, with a crossover at tkjJk ~ /i^sT, with 
width in huk/ksT equal to that of the undamped harmonic lattice. 

The width of the crossover region is defined in terms of the ratio fiWk/kBT (see Fig. |I|). 
The phonons excluded from the finite system have linear dispersion with an upper bound 
of UJk = skc = 27rsp^^^/N^^^, where s is the speed of sound. By choosing typical values for 
temperature, density and sound speed, we can estimate the width of the crossover region as 
a function of the number of particles N. For example, if we assume the speed of sound is 



s = 4.0 X 10^^ A/sec, we find that for temperature T = 20.0 K and density p = 0.055 A" 
-^crossover ~ 10^ — 10^. A simple power law fit to system sizes in this range is inaccurate. 
If we assume a harmonic spectrum, we can use Eq. ([16|) to write 



6{u^) oc dk P A{k) oc dk fccoth( 






'l-exp{-BN-^/^ 

BN-y^ 




0{N-'), (19) 

as 




where B = 2TcPhsp^^^. The functional form of (m^) can now be written 

l + l°g( |)v-V3 ' jj^ (20) 

The parameter B is not effectively constrained by the PIMC values of (m^) versus A^, and 
therefore cannot be accurately determined by fitting. Instead, by using physically reasonable 
values of 5, we can use least-squares fitting (of A and (m^)oo) to extrapolate our PIMC data 
to the thermodynamic limit and get values for (m^)oo (see Table |I[) The speed of sound s 
can be obtained from experimental measurements of isothermal compressibility Kitll, using 
the relationship 

s' = ^, (21) 

prriKT 

However, we obtained poor fits using the value of B obtained in this manner. The lack of 
self-consistency with physical parameters shows that the undamped harmonic lattice used to 
derive Eq. ( pil| ) is insufficient for describing the mean-squared displacement in solid helium. 
A larger range of B was required for reasonable fitting, both due to anisotropy in s and the 
approximations used in calculating the functional form Eq. (|20|) . We were able to obtain 
reasonable fitting by using values of B which correspond to s = 4000 ± 2000 m/sec, which 
is larger than experimental values by about a factor of 2 on average. For systems away from 
melting, the fitting was insensitive to the value of B, and the propagated error from B was 
on the same order of magnitude as the fitting error, in estimating (m^)oo- The reduced x^ 
for systems away from melting was typically between 0.1 and 1.0. For systems near melting, 
the fitting was much more sensitive to the estimated range of B, particularly the lower 
bound, and the errorbar on (m^)oo was completely determined by the uncertainty in B. The 
reduced x^ &!■ these systems was usually between 5 and 10. A more accurate functional 
form, perhaps one which takes the phonon linewidth into account, would greatly improve 
the accuracy of extrapolating to the thermodynamic limit for these systems. 

Our PIMC calculations agree with all direct scattering measurements of (m^), when 
extrapolated to the thermodynamic limit. We are able to confirm both computational 
and experimental methods to an accuracy of 5% in the mean-squared displacement. Near 
melting, the accuracy of the PIMC values was considerably reduced, due to uncertainty in the 
functional form of the finite-size effects. PIMC values were generally higher than the indirect 
scattering measurement si^. Because the indirect measurements assumed contributions from 
single-phonon processes only, this discrepancy gives evidence for the importance of multi- 
phonon processes in solid helium. 
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IV. TEMPERATURE DEPENDENCE 



Using Eq. (|r^), it is straightforward to show that the harmonic approximation predicts 
(m^)oo to have temperature dependence of the form 

{u')oo = {uYoo=' + CT'. (22) 

However, our extrapolated values of (m^) appear to fit to a T^ power law, rather than T^ 
for 5 K < T < 20 K (see Fig. [^]). With less than a decade in temperature, it is difficult 
to draw any direct conclusions from this, although it may indicate that harmonic theory is 
insufficient to accurately describe solid helium, and that higher order terms dominate the 
temperature dependence. 



V. NON-GAUSSIAN CORRECTIONS TO (u^) 

We have determined the deviation of the density from a Gaussian distribution by two 
methods. The first was fitting ln{S{k)) to a polynomial in k"^. As shown in Eq. (|^), the 
linear term is (u^) , the quadratic term is proportional to the kurtosis n. We also directly 
calculated the kurtosis in the (111) and (100) directions, using Eq. (|). The kurtosis was 
found to be non-zero and anisotropic in the fee solid hehum systems we studied. 

Shown are graphs with kurtosis in the (100) direction as a function of density and tem- 
perature. The kurtosis is roughly twice as large in systems near melting. Away from melting, 
the kurtosis appears to be independent of both temperature and density. 

We found that k is independent of A^, for values of A^ > 500. For smaller system 
sizes, the finite-size effects are large, but drop off quickly with increasing N. However, at 
the wavevectors k < 9 A^^, the effect of the kurtosis is only a few percent, making this 
term difficult to observe in scattering experimentsiSlMEjy'y. Vitiello et allis have computed 
kurtoses of 0.051 and 0.042 at zero temperature using Shadow wave functions, for molar 
volumes of 20.5 and 18.3. These values are consistent with the values computed here. 

VI. CONCLUSIONS 

PIMC simulations of the Debye- Waller factor in solid helium agree with experimental 
results to better than 5% accuracy, indicating that the assumed potential, the computational 
methods, and the experimental analysis are correct within the stated errors. We determined 
the first non-Gaussian contribution, a directionally-dependent kurtosis. The finite-size ef- 
fects were found to be in a crossover region between the classical and quantum scaling limits, 
and hence a power-law dependence was insufficient for extrapolation to the thermodynamic 
limit. The harmonic approximation gave an approximate functional form for (m^) which was 
used to extrapolate finite values to the thermodynamic limit. The extrapolated values agree 
with all available direct scattering measurements, although the extrapolation became more 
sensitive with increased proximity to the melting transition, with correspondingly larger un- 
certainties. Higher extrapolation accuracy could be obtained with a more accurate scaling 
form. The effective temperature dependence of (m^) appears to be closer to T^, rather than 
the T^ predicted by harmonic theory. 
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TABLES 



TABLE L PIMC results, (u )oo was estimated by fitting finite PIMC data to the function 
given in Eq. (P0|), assuming a sound speed s = 4000 =b 2000 m/see. The listed uneertainty in 
{u^)oo represents the range of fitted values corresponding to the uneertainty in s. The errors in the 
energies represent the statistical uncertainties in the final digit, while the experimental errors are 
the total uncertainties as estimated by the authors, kiqo represents the estimated kurtosis in the 
(100) direction, k was estimated for the largest system size available, usually A^ = 1372. 



type 


Kn(cm3) 


T(K) 


Texpt (K) 


(^2)^(10-2 A2) 


(n2),,p,(10-2 A2) 


/^lOO 


-C/kinotic (K 


) E,,, (K) 


fee ^He 


10.98 


20.00 


20.25 


9.77(39) 


9.99(27)^ 


0.09 


80.30(2) 


50.24(2) 


fee ^He 


10.98 


17.78 




9.15(32) 




0.08 


79.10(4; 


48.02(5) 


fee ^He 


10.98 


16.84 




8.71(23) 




0.07 


79.15(6) 


47.52(6) 


fee ^He 


10.98 


16.00 




8.76(21) 




0.07 


78.51(4) 


46.64(4) 


fee ^He 


10.98 


15.24 




8.49(18) 




0.06 


78.61(6) 


46.44(6) 


fee ^He 


10.98 


13.33 




8.36(16) 




0.06 


77.63(4) 


45.08(4) 


fee ^He 


10.98 


10.00 




8.02(8) 




0.06 


77.06(2) 


43.99(2) 


fee ^He 


10.98 


8.00 




7.89(7) 




0.06 


76.87(4) 


43.67(3) 


fee %e 


10.98 


5.00 




7.79(5) 




0.06 


76.80(2) 


43.50(2) 


fee ^He 


10.39 


24.60 


24.40 


8.89(53) 


8.23(40)^ 


0.09 


90.41(4) 


68.71(5) 


fee ^He 


10.02 


26.67 


25.94 


8.29(61) 


8.43(18)^ 


0.09 


96.48(4) 


81.53(4) 


fee ^He 


9.02 


35.56 


38.00 


6.91(94) 


5.66= 


0.09 


118.14(7) 


133.74(8) 


fee ^He 


9.02 


22.86 




5.77(21) 




0.04 


110.68(5) 


120.77(5) 


fee ^He 


9.43 


21.33 




6.31(21) 




0.05 


109.14(6) 


114.16(6) 


fee ^He 


9.97 


26.67 


28.00 


8.14(58) 


6.93= 


0.08 


101.33(5) 


94.26(5) 


fee ^He 


9.97 


18.82 


19.00 


7.04(20) 


6.20= 


0.05 


98.82(6) 


88.03(6) 


fee ^He 


11.54 


17.78 


18.13 


11.41(40) 


11.43(11)^ 


0.10 


86.14(3) 


58.53(3) 


fee ^He 


11.54 


10.00 




10.00(12) 




0.08 


84.15(4) 


54.68(4) 


fee ^He 


11.54 


5.00 




9.72(5) 




0.12 


83.88(5) 


54.22(4) 


fee 3He 


10.98 


17.78 




9.87(28) 




0.09 


93.18(5) 


70.39(5) 


fee 3He 


10.00 


29.09 




9.07(76) 




0.09 


113.74(5) 


110.00(6) 


hep ^He 


12.12 


14.55 


14.23 


11.17(13) 


11.25(28)^ 




66.44(4) 


27.95(4) 


hep ^He 


12.12 


11.85 


12.00 


10.24(13) 


10.26(17f 




65.56(6) 


26.16(6) 


hep ^He 


12.12 


5.00 




9.52(5) 






64.37(3) 


24.15(3) 


hep ^He 


15.72 


5.71 


5.80 


17.40(13)<= 


17.31(18)'^ 




40.66(3) 


0.78(3) 


hep ^He 


10.98 


5.00 




7.78(3) 






76.74(4) 


43.47(3) 


hep ^He 


11.90 


16.84 


16.81 


11.84(24) 


11.96(26)^ 




82.17(6) 


52.13(6) 


hep ^He 


12.81 


12.31 


12.54 


13.26(19) 


13.43(27f 




71.14(4) 


36.37(4) 



'^Arms and Simmonst3. Direet x-ray measurements. 
Venkataraman and Simmonsll3. Direet x-ray measurements. 
=Thomlinson, Eekert and ShiraneEJ. Indirect neutron measurements . 
<^Stass.s, Kh^tamian, and KhnJi. Direet neutron measurements. 

'^ExtraDolated using speed of sound s = 1000 it 200 m/sec, based on experimental data near this 
densityii^. 
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FIGURES 

FIG. 1. Numerical calculation of A{k) as a function of hwh/ksT, for a damped harmonic 
oscillator with linewidth 7^ = 0.2a;fc. Shown is the crossover between the quantum limit (dashed 
line) and the classical limit (dotted line), which allows one to estimate the range of system sizes 
affected by the crossover region. 

FIG. 2. Extrapolation of finite PIMC simulation data to the thermodynamic limit, for fee 
^He , Vm = 10.98 cm^, T=20.0 K. {v?) was calculated by direct averaging (open circles) and by 
fitting to S{k) data (open triangles). Both methods always agreed within statistical error bars. 
Least-squares fitting was used to fit the directly averaged data to Eq. (|20|). The experimental data 
point (solid circle) is from Arms and Simmons. 

FIG. 3. Temperature dependence of {v?)oo , using the extrapolated values from Table |. Shown 
are values for fee ^He, Vm = 10.98 cm^ (solid circles), and fee ^He, Vm = 11.54 cm^ (solid squares). 

FIG. 4. ^ln(^^^) versus fc^ for fee ^He , Vm = 10.98 cm^, N=864, T=20.0 K. The kurtosis 
K is given by the slope, and agrees with direct calculations of {u^) / {u^)"^ — 3. In the (100) direction 
(dashed line), the fitted kurtosis was 0.11±0.03, while the direct value was 0.09±0.03. In the (111) 
direction (dotted line), the fitted kurtosis was 0.05 it 0.03, while the direct value was 0.02 ± 0.03. 

FIG. 5. Kurtosis k, vs. molar volume Vm and temperature T, in the (100) direction, for fee '^He. 
The kurtosis is noticeably larger near the experimental melting line, but is otherwise independent 
of temperature and density. 
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